// Diffractor.cpp: implementation of the Diffractor class.
//
//////////////////////////////////////////////////////////////////////

#define PI   3.1415926535897932385

#include "Diffractor.h"

#include <fstream>  // used to get ifstream
#include <sstream>
#include <algorithm>

extern "C"
{
#include "../Atominfo/atominfo.h"
}

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////


Diffractor::Diffractor(char *DiffractorFile,
                                           Position *W1, Position *W2, Position *W3, char **atomLabel, int nAtom)
{

  NAtom = nAtom;  // note NAtom is redefined in Diffractor_read_old_style()

  string filetype(DiffractorFile);

  size_t found = filetype.find(".hcv_dash");
  if (found != string::npos)
    Diffractor::Diffractor_read_hcv(DiffractorFile, W1, W2, W3, atomLabel);
  else
    Diffractor::Diffractor_read_old_style(DiffractorFile, W1, W2, W3, atomLabel);
}




Diffractor::~Diffractor()
{
delete  Data;

deleteA(Intensity);
deleteA(ReAmplitude);
deleteA(ImAmplitude);
deleteA(H);
deleteA(K);
deleteA(L);
deleteA(DH);
deleteA(DK);
deleteA(DL);
deleteM(FiniteFactor);
deleteM(FDeriv);

deleteSymmA(SinX, MaxH);
deleteSymmA(CosX, MaxH);
deleteSymmA(SinY, MaxK);
deleteSymmA(CosY, MaxK);
deleteSymmA(SinZ, MaxL);
deleteSymmA(CosZ, MaxL);
}


//////////////////////////////////////////////////////////////////////
// functions
//////////////////////////////////////////////////////////////////////

double *Diffractor::Build(Position *AtomPosition)
{
FourierTransform(AtomPosition);

for(int n=1;n<=NIntensity;n++)
        Intensity[n] = ReAmplitude[n]*ReAmplitude[n] + ImAmplitude[n]*ImAmplitude[n];

return (Intensity);
}


////////////////////////////////////////
// private methods
////////////////////////////////////////

void Diffractor::Diffractor_read_hcv(char *DiffractorFile,
                                           Position *W1, Position *W2, Position *W3, char **atomLabel)
{

  // See comment top for why this code commented out

  //cout << "\nDiffractor_read_hcv method code commented out for now\n";
	cout << "\nReading from hcv_dash format\n";
  //exit(1);




  // Number of off-diagonal weight matrix elements read in for each line
  // of the 'new' style .hcv file.
  // Notice when a multiplet SPLIT it into singlets
  vector<int> num_cor_element;



  ///* M-EDIT
  // read diffraction data into attribute data which has
  // the structure:

  int*    h;                   // pointer to an array of h's
  int*    k;                     
  int*    l;
  double* ave_intensity;       // pointer to an array of the likelihood expectation values of the intensities
  double **inv_cov;            // pointer to inverse of the likelihood covariance matrix of the intensities
  int N;   // the number of intensities
  //*/



  ///* M-EDIT
  populate_store_data(DiffractorFile, num_cor_element);


  FracX = (2.0*PI)*(*W1);
  FracY = (2.0*PI)*(*W2);
  FracZ = (2.0*PI)*(*W3);

  M = TPMatrix(FracX,FracY,FracZ);

  NIntensity = data.N;

  Data = new DataHolder(NIntensity);

  newA(Intensity,               NIntensity);
  newA(H,                               NIntensity);
  newA(K,                               NIntensity);
  newA(L,                               NIntensity);
  newA(DH,                      NIntensity);
  newA(DK,                      NIntensity);
  newA(DL,                      NIntensity);
  newA(ReAmplitude,     NIntensity);
  newA(ImAmplitude,     NIntensity);
  newM(FiniteFactor,    NAtom, NIntensity);
  newM(FDeriv,          NAtom, 3);

  // to calculate average form factor decay

  vector<double> averageDecay;


  double weightMatrixOffDiagonalCutoff = 20.0;
  
  // populate Data.Hessian from inv_cov array
  // as read in from hcv file, where diagonal elements
  // are i/sigma and off diagonal elements percentages

  for (int i = 1; i <= NIntensity; i++)
  {
    Data->Hessian[i][i] = data.inv_cov[i][1]*data.inv_cov[i][1];
  }

  // filling in the off diagonal elements

  for (int i = 1; i <= NIntensity; i++)
  {
    for (int j = i+1; j-i+1 <= num_cor_element[i]+1; j++)
    {  
      if (data.inv_cov[i][j-i+1] > weightMatrixOffDiagonalCutoff)
      {
          // convert percentage inv_cov into absolute hessian
          Data->Hessian[i][j] = data.inv_cov[i][j-i+1]
            * sqrt(Data->Hessian[i][i] * Data->Hessian[j][j]) / 100.0;

          Data->Hessian[j][i] = Data->Hessian[i][j];
      }
    }
  }



  MaxH = 0;
  MaxK = 0;
  MaxL = 0;

  for(int m=1;m<=NIntensity;m++)
        {
          Data->Average[m] = data.ave_intensity[m];
        
          H[m] = data.h[m];
          K[m] = data.k[m];
          L[m] = data.l[m]; 


          DH[m] = static_cast<double>(H[m]);
          DK[m] = static_cast<double>(K[m]);
          DL[m] = static_cast<double>(L[m]);



          if(abs(H[m])>MaxH) MaxH = abs(H[m]);
          if(abs(K[m])>MaxK) MaxK = abs(K[m]);
          if(abs(L[m])>MaxL) MaxL = abs(L[m]);
        }

  newSymmA(SinX, MaxH);
  newSymmA(CosX, MaxH);
  newSymmA(SinY, MaxK);
  newSymmA(CosY, MaxK);
  newSymmA(SinZ, MaxL);
  newSymmA(CosZ, MaxL);


  SinX[0] = 0.0;                CosX[0] = 1.0;
  SinY[0] = 0.0;                CosY[0] = 1.0;
  SinZ[0] = 0.0;                CosZ[0] = 1.0;


  // creating your own FiniteFactors

  for(int m=1;m<=NIntensity;m++)
  {
     Position temp_pos = (*W1)*H[m]+(*W2)*K[m]+(*W3)*L[m];
     double stol2 = temp_pos*temp_pos / 4.0;

           for(int n=1;n<=NAtom;n++)
     {  
        const T_SF_WK95_CAA *AtomScatFac;
        AtomScatFac = FindSF_WK95_CAA(atomLabel[n], 0);
        FiniteFactor[n][m] = CalcSF_WK95_CAA(AtomScatFac, stol2);
     }

     // optionally include a debye-waller fall off if this has not 
     // already been done in hcv file. 
     // In general depending on data this may have little effect
     // The number 2.42 seems a bit arbitray and of course will
     // be heavily influenced by temperature
           for(int n=1;n<=NAtom;n++)
     {  
       FiniteFactor[n][m] *= exp(-2.42*stol2);
     }
  }


//*/
}


void Diffractor::Diffractor_read_old_style(char *DiffractorFile,
                                           Position *W1, Position *W2, Position *W3, char **atomLabel)
{
  int    i, j, n, m, NElement, NSpecial;
  double d1, d2, d3, dummy;

  SAFEifstream infile(DiffractorFile);


  FracX = (2.0*PI)*(*W1);
  FracY = (2.0*PI)*(*W2);
  FracZ = (2.0*PI)*(*W3);

  M = TPMatrix(FracX,FracY,FracZ);

  infile >> NIntensity >> NElement >> NAtom >> NSpecial;

  Data = new DataHolder(NIntensity);

  newA(Intensity,               NIntensity);
  newA(H,                               NIntensity);
  newA(K,                               NIntensity);
  newA(L,                               NIntensity);
  newA(DH,                      NIntensity);
  newA(DK,                      NIntensity);
  newA(DL,                      NIntensity);
  newA(ReAmplitude,     NIntensity);
  newA(ImAmplitude,     NIntensity);
  newM(FiniteFactor,    NAtom, NIntensity);
  newM(FDeriv,          NAtom, 3);


  for(m=1;m<=NElement;m++)
          {
          infile >> i >> j >> d1 >> d2 >> d3;
        
          if(i>NIntensity || j>NIntensity)
                  {
                  cerr << "Problem with Hessian elements in *.inf!\n\n";
                  exit(1);
                  }

          if(i==j)      
                  Data->Hessian[i][i] = d1;
          else
                  {
                  Data->Hessian[i][j] = d1/2.0;
                  Data->Hessian[j][i] = d1/2.0;
                  }
          }



  MaxH = 0;
  MaxK = 0;
  MaxL = 0;

  for(m=1;m<=NIntensity;m++)
          {
          infile >> d1 >> d2 >> d3 >> Data->Average[m];
        
          //for reading in special .inf
  /*    int int_d1, int_d2, int_d3;
          infile >> i >> j >> int_d1 >> int_d2 >> int_d3;


          H[m] = int_d1;
          K[m] = int_d2;
          L[m] = int_d3;

          DH[m] = double(H[m]);
          DK[m] = double(K[m]);
          DL[m] = double(L[m]);*/

          DH[m] = ceil(d1/(2.0*PI) - 0.5);
          DK[m] = ceil(d2/(2.0*PI) - 0.5);
          DL[m] = ceil(d3/(2.0*PI) - 0.5);

          H[m] = int(DH[m]);
          K[m] = int(DK[m]);
          L[m] = int(DL[m]); 

          if(abs(H[m])>MaxH) MaxH = abs(H[m]);
          if(abs(K[m])>MaxK) MaxK = abs(K[m]);
          if(abs(L[m])>MaxL) MaxL = abs(L[m]);
          }


  newSymmA(SinX, MaxH);
  newSymmA(CosX, MaxH);
  newSymmA(SinY, MaxK);
  newSymmA(CosY, MaxK);
  newSymmA(SinZ, MaxL);
  newSymmA(CosZ, MaxL);


  SinX[0] = 0.0;                CosX[0] = 1.0;
  SinY[0] = 0.0;                CosY[0] = 1.0;
  SinZ[0] = 0.0;                CosZ[0] = 1.0;


  for(m=1;m<=NIntensity;m++)
          for(n=1;n<=NSpecial;n++)
                  infile >> dummy;

  // for reading in special .inf

  /*for(m=1;m<=NIntensity;m++)
          for(n=1;n<=NSpecial;n++)
                  infile >> fisse;*/

  for(m=1;m<=NIntensity;m++)
          for(n=1;n<=NAtom;n++)
                  infile >> FiniteFactor[n][m];


}




// ***********************************************************************
// This function checks if the file infilename has the format
//              h(int) k(int) l(int) Int(double) Sigma(double) Counter(long) ...
//              CorToNext(double) CorToSecond(double) etc.      
// And returns the number of lines and the number of doubles in infilename through
// its first two arguments
// ***********************************************************************

void Diffractor::get_information_infile(long *number_of_lines_ptr, long *num_extra_lines_ptr, 
                                                        long *highest_multiplet_ptr, char infilename[])
{
        long highest_multiplet_dummy;   // used in determining highest_multiplet

        FILE *infile_fp;                // where the file infilename is kept
        
  int dummy_h, dummy_k, dummy_l;                // used to read in numbers from the
        double dummy_int, dummy_sigma;          // file infilename
        long dummy_counter;                                     //      
                
        int previous_counter;                           // used to check for doubles
        
        char line[5000];                                                // to read in first line of infilename
        
        
        // open infilename
        
        infile_fp = Fopen(infilename,"r");

        
        // read in first line   
        
        if (fgets(line, 5000, infile_fp) == NULL) 
        {
                printf("\nError reading first line of %s\n", infilename);

                exit(1);
        }



        // checks that the format of the input file is not completely crazy

        if ( 6 == sscanf(line, "%d%d%d%lf%lf%ld", &dummy_h, &dummy_k, &dummy_l, 
                &dummy_int, &dummy_sigma, &dummy_counter) )
        {

        }
        else
        {
                printf("The input format of %s is invalid\n", infilename);
                printf("Please see manual.doc explaining the required format of the input file\n");

                exit(1);
        }
        
        
        // move the file pointer associated with infile_fp back to the beginning of the file
        
        fseek(infile_fp, 0L, SEEK_SET);
        
        
        // initiate previous_counter to dummy number
        
        previous_counter = -1;
        
        
        // initiate number of lines and multiplet
        
        *number_of_lines_ptr = 0;
        *num_extra_lines_ptr = 0;
        *highest_multiplet_ptr = 0;
        highest_multiplet_dummy = 0;
        
        // find the number of lines etc
        
        while ( fgets(line, 5000, infile_fp) != NULL )
        {
                if (count_number_in_line(line) == 0)
                {
                        break;
                }
        
                // increase number of lines by one
                
                ++(*number_of_lines_ptr);
                
                
                // check for the presence of a multiplet
                
                sscanf(line, "%d%d%d%lf%lf%ld", &dummy_h, &dummy_k, &dummy_l, 
                        &dummy_int, &dummy_sigma, &dummy_counter);
                
                if (previous_counter == dummy_counter)
                {
                        // increase number of doublets by one
                        
                        ++(highest_multiplet_dummy);
                        ++(*num_extra_lines_ptr);
                }
                else
                {
                        if (highest_multiplet_dummy > *highest_multiplet_ptr)
                                *highest_multiplet_ptr = highest_multiplet_dummy;

                        // put highest_multiplet_dummy back to zero again

                        highest_multiplet_dummy = 0;

                        // set previous_counter to current counter
                
                        previous_counter = dummy_counter;
                }
        }
    

        // close file
        
        fclose(infile_fp);
}


FILE* Diffractor::Fopen(char *fname, char *fmode) {
  FILE *fp;
  
  if ((fp=fopen(fname, fmode)) == NULL)
  {
    cout << "Couldn't open file " << fname << " in mode " << fmode << endl;
  }
  return(fp);
}


///////////////////////////////////
// Read in a 'new' type Dash hcv file.
//
// However split any multiplet into singlet. 
// Assume multiplet have the format:
//     h k l ..... 32 100 4
//     h k l ..... 32 4
///////////////////////////////////
void Diffractor::populate_store_data(char* DiffractorFile, vector<int>& num_cor_element)
{               
        int counter;                    // the counter in the data file

        long num_extra_lines;                           // number of extra lines from multiplets

        long highest_multiplet;    // highest multiplet in input data file
                                                                                                                 // highest_multiplet 0 = only singlets
                                                                                                                 // highest_multiplet 1 doublet and so on

        long number_of_lines;                   // equal to num_refined_reflections + num_extra_lines

  long PAWLEY_COR_ELEMENTS = 20;


        // checking and getting number of refined reflections etc from infile
        
        get_information_infile(&number_of_lines, &num_extra_lines, &highest_multiplet, 
                DiffractorFile);
                

        // open infile which will be read below 
        
        ifstream data_file(DiffractorFile); 
        
        
        data.N = number_of_lines;   

  num_cor_element.resize(data.N+1); // +1 to conform to Numerical recipes standard

        newM(data.inv_cov, data.N, PAWLEY_COR_ELEMENTS);

  {
  for (int i = 1; i <= data.N; i++)
    for (int j = 1; j <= PAWLEY_COR_ELEMENTS; j++)
      data.inv_cov[i][j] = 0.0; // done since the code calling this method make use of none populated elements are zero
  } // somehow microsoft visual studio gets confused here and don't properly clear i and j

  newA(data.ave_intensity, data.N);
        
        newA(data.h, data.N);
        newA(data.k, data.N);
        newA(data.l, data.N);


        // initialize last_counter to dummy number


        int i = 0;   // pointing at first peak in infilename

  string aLine;

        while( getline(data_file, aLine, '\n') )
  {
    int number = count_number_in_line(aLine.c_str());

    if (number == 0)
      break;

    i++;

    num_cor_element[i] = number - 6; // h,k,l,I,sigma,count not counted


    stringstream ins(aLine);

                // read in h,k,l
                
                ins >> data.h[i] >> data.k[i] >> data.l[i];
                
                ins >> data.ave_intensity[i] >> data.inv_cov[i][1];


    // read in the data counter
                
                ins >> counter;
                
                
                        
                for(int j=2;j<=num_cor_element[i]+1;j++)
    {
                        ins >> data.inv_cov[i][j];

      // In principle the input file should not contain -100 or 100 correlation 
      // elements - since if two peaks are 100% correlated than better to treat
      // these as a one refineable. Further to avoid rule out the possible that
      // a matrix inversion struggle then do the stuff below

      if (data.inv_cov[i][j] < -100.0001 || data.inv_cov[i][j] > 100.0001)
      {
        printf("\n\n***********************************************************\n\n");
                          printf("ERROR IN INPUT FILE:\n");     
                          printf("  Correlation element for reflection %d %d %d\n", 
          data.h[i], data.k[i], data.l[i]);
                          printf("  is larger than 100 percent, which is clearly nonsense\n.");
                          printf("***********************************************************\n\n");

                          exit(1);
      }
    }
  }

        // double check that infilename has been stored correctly in data[0]
        
        if (i != data.N)
        {
                printf("Error in reading %s\n",
                        DiffractorFile);
                exit(1);
        }

}


// return number of items in line
int Diffractor::count_number_in_line(const char* line)
{
                int i = 0; // plays the role of stepping through the string line
                int j = 0; // plays the role of counting the number of numbers in line
                do
                {
                        i++;
                        if ((line[i] == ' ' || line[i] == '\t') && (line[i-1] != ' ' && line[i-1] != '\t'))
                        {
                                j++;
                        }
                        if (line[i] == '\0') break;
                } while (line[i] != '\n');
                if (line[i-1] != ' ' && line[i-1] != '\t') j++; 
                
    return j;
}


